Data Read In

CSLAP<-read.csv("CSLAP_Dataset_09232019.csv", header=TRUE)
CSLAP$Sample_Year<-as.factor(CSLAP$Sample_Year)
CSLAP$Sample_Month<-as.factor(CSLAP$Sample_Month)

#Split by `Info_Type`
OWCSLAP<-CSLAP[CSLAP$Info_Type == "OW",]
BSCSLAP<-CSLAP[CSLAP$Info_Type == "BS",]
SBCSLAP<-CSLAP[CSLAP$Info_Type == "SB",]
noSBCSLAP<-CSLAP[CSLAP$Info_Type != "SB",]

redCSLAP<-read.csv("redCSLAP.csv")

redOWCSLAP<-redCSLAP[redCSLAP$Info_Type == "OW",]
redBSCSLAP<-redCSLAP[redCSLAP$Info_Type == "BS",]
redSBCSLAP<-redCSLAP[redCSLAP$Info_Type == "SB",]
rednoSBCSLAP<-redCSLAP[redCSLAP$Info_Type != "SB",]

Mixed Effects Models for the Global Dataset (n=69 lakes)

Open Water Total Phosphorus

getwd()
## [1] "/Users/victoriafield/Library/Mobile Documents/com~apple~CloudDocs/Masters Thesis/Chapter Two CSLAP Manuscript/R - data and outputs"
#Select out data of interest
df <- OWCSLAP %>% dplyr::select(TP_mg.L, Dreissenids, CA.SA, Mean_Depth_m, Sample_Year, Sample_Month, LakeID)
df <- df[complete.cases(df), ]

#Check distribution 
##gamma
gamma <- fitdistr(df$TP_mg.L+0.01, "gamma")
qqp(df$TP_mg.L+0.01, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])

## [1] 2665 2659
##log-normal 
qqp(df$TP_mg.L, "lnorm")

## [1] 2665 2659

Model 1 - Full Model, TP log transformed

fit1<-lmer(log(TP_mg.L+.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1|Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df)

summary(fit1)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(TP_mg.L + 0.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m +  
##     (1 | Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 |  
##     Sample_Year:LakeID)
##    Data: df
## 
## REML criterion at convergence: -756.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6502 -0.4731 -0.0937  0.3083  9.7084 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 0.0175696 0.13255 
##  LakeID             (Intercept) 0.0256021 0.16001 
##  Sample_Month       (Intercept) 0.0001905 0.01380 
##  Sample_Year        (Intercept) 0.0022388 0.04732 
##  Residual                       0.0338332 0.18394 
## Number of obs: 2702, groups:  
## Sample_Year:LakeID, 349; LakeID, 69; Sample_Month, 6; Sample_Year, 6
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -3.721277   0.082613  86.712111 -45.045  < 2e-16 ***
## DreissenidsUninvaded  -0.049672   0.059274 111.130287  -0.838 0.403825    
## log(CA.SA)             0.018196   0.021658  65.195196   0.840 0.403900    
## Mean_Depth_m          -0.022151   0.005685  66.781005  -3.896 0.000229 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.678              
## log(CA.SA)  -0.545  0.058       
## Mean_Dpth_m -0.294  0.014 -0.206
plot(fit1,col=df$LakeID) #plots model residuals grouped by lake

plot(fit1,col=df$LakeID, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

plot(fit1,col=df$Sample_Year) #plots model residuals grouped by year 

plot(fit1,col=df$Sample_Year, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

qqnorm(resid(fit1))
qqline(resid(fit1))

From observation numbers, we can see that those points on the far right of each plot correspond to the following lakes and years: Guilford Lake, East Caroga, Little Long Pond, and Galway Lake. And Year 2015 (blue), Year 2012 (black), Year 2017 (Pink), and Year 2016 (Cyan)

Model 2 - removing “troublesome” lakes

#Select out data of interest
df2 <- df[df$LakeID != "0601GUI0188",]
df2 <- df2[df2$LakeID != "1201ECA0697", ]
df2 <- df2[df2$LakeID != "1701LLO0708", ]
df2 <- df2[df2$LakeID != "1201GAL0563", ]

fit2<-lmer(log(TP_mg.L+.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1|Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df2)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0170704 (tol = 0.002, component 1)
summary(fit2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(TP_mg.L + 0.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m +  
##     (1 | Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 |  
##     Sample_Year:LakeID)
##    Data: df2
## 
## REML criterion at convergence: -1234.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0533 -0.4788 -0.1202  0.3204  7.7073 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 0.0070795 0.08414 
##  LakeID             (Intercept) 0.0232003 0.15232 
##  Sample_Month       (Intercept) 0.0002114 0.01454 
##  Sample_Year        (Intercept) 0.0011476 0.03388 
##  Residual                       0.0289562 0.17017 
## Number of obs: 2515, groups:  
## Sample_Year:LakeID, 325; LakeID, 65; Sample_Month, 6; Sample_Year, 6
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -3.735775   0.073489  88.292693 -50.834  < 2e-16 ***
## DreissenidsUninvaded  -0.053305   0.050551 136.143312  -1.054 0.293530    
## log(CA.SA)             0.017758   0.020410  61.089340   0.870 0.387659    
## Mean_Depth_m          -0.021598   0.005322  62.033773  -4.058 0.000141 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.644              
## log(CA.SA)  -0.571  0.054       
## Mean_Dpth_m -0.304  0.006 -0.213
## convergence code: 0
## Model failed to converge with max|grad| = 0.0170704 (tol = 0.002, component 1)
plot(fit2,col=df2$LakeID) #plots model residuals grouped by lake

plot(fit2,col=df2$LakeID, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

plot(fit2,col=df2$Sample_Year) #plots model residuals grouped by year 

plot(fit2,col=df2$Sample_Year, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

qqnorm(resid(fit2))
qqline(resid(fit2))

Removing lakes improves residual plot, but not QQ plot

Model 3 - removing “troublesome” Year 2015

#Select out data of interest
df3 <- df[df$Sample_Year != "2015",]


fit3<-lmer(log(TP_mg.L+.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1|Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df3)

summary(fit3)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(TP_mg.L + 0.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m +  
##     (1 | Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 |  
##     Sample_Year:LakeID)
##    Data: df3
## 
## REML criterion at convergence: -1128.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0512 -0.4997 -0.1188  0.3060 10.4227 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 5.260e-03 0.072526
##  LakeID             (Intercept) 2.742e-02 0.165578
##  Sample_Month       (Intercept) 2.429e-05 0.004929
##  Sample_Year        (Intercept) 1.690e-03 0.041110
##  Residual                       2.903e-02 0.170379
## Number of obs: 2272, groups:  
## Sample_Year:LakeID, 292; LakeID, 69; Sample_Month, 6; Sample_Year, 5
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -3.764171   0.077047  96.364793 -48.855  < 2e-16 ***
## DreissenidsUninvaded  -0.049569   0.051430 169.912048  -0.964 0.336509    
## log(CA.SA)             0.023500   0.021539  65.718248   1.091 0.279229    
## Mean_Depth_m          -0.019589   0.005633  66.493563  -3.477 0.000898 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.630              
## log(CA.SA)  -0.571  0.048       
## Mean_Dpth_m -0.311  0.010 -0.204
plot(fit3,col=df3$LakeID) #plots model residuals grouped by lake

plot(fit3,col=df3$LakeID, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

plot(fit3,col=df3$Sample_Year) #plots model residuals grouped by year 

plot(fit3,col=df3$Sample_Year, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

qqnorm(resid(fit3))
qqline(resid(fit3))

Removing 2015 improves neither residuals nor QQ

Model 3 - removing “troublesome” lakes AND 2015

#Select out data of interest
df4 <- df[df$LakeID != "0601GUI0188",]
df4 <- df4[df4$LakeID != "1201ECA0697", ]
df4 <- df4[df4$LakeID != "1701LLO0708", ]
df4 <- df4[df4$LakeID != "1201GAL0563", ]
df4 <- df4[df4$Sample_Year != "2015",]

fit4<-lmer(log(TP_mg.L+.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1|Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df4)

summary(fit4)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(TP_mg.L + 0.01) ~ Dreissenids + log(CA.SA) + Mean_Depth_m +  
##     (1 | Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 |  
##     Sample_Year:LakeID)
##    Data: df4
## 
## REML criterion at convergence: -1376.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.4509 -0.5243 -0.1206  0.3339  8.2501 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 3.534e-03 0.059448
##  LakeID             (Intercept) 2.470e-02 0.157154
##  Sample_Month       (Intercept) 1.454e-05 0.003813
##  Sample_Year        (Intercept) 1.275e-03 0.035706
##  Residual                       2.515e-02 0.158588
## Number of obs: 2114, groups:  
## Sample_Year:LakeID, 272; LakeID, 65; Sample_Month, 6; Sample_Year, 5
## 
## Fixed effects:
##                        Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)           -3.767159   0.072209  91.587601 -52.170  < 2e-16 ***
## DreissenidsUninvaded  -0.048822   0.046906 174.969528  -1.041 0.299381    
## log(CA.SA)             0.017846   0.020804  61.360106   0.858 0.394355    
## Mean_Depth_m          -0.018813   0.005418  61.964622  -3.472 0.000947 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.606              
## log(CA.SA)  -0.584  0.046       
## Mean_Dpth_m -0.315  0.004 -0.213
plot(fit4,col=df4$LakeID) #plots model residuals grouped by lake

plot(fit4,col=df4$LakeID, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

plot(fit4,col=df4$Sample_Year) #plots model residuals grouped by year 

plot(fit4,col=df4$Sample_Year, id=0.5, idLabels=~.obs, cex=.5) #plots model residuals with labels (observation number) for values outside Normal confidence limits 

qqnorm(resid(fit4))
qqline(resid(fit4))

Removing “bad” lakes and Year 2015 improves residual plot the most, but the QQ line is still wonky

Alternate distributions using Generalized Linear Mixed Effects Model

Model 1 - Gaussian Family with Log Link

#model transformed
Gfit<-glmer(TP_mg.L+.01 ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1 | Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df, family=gaussian(link="log"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.433841 (tol = 0.001, component 1)
summary(Gfit)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: gaussian  ( log )
## Formula: TP_mg.L + 0.01 ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1 |  
##     Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 | Sample_Year:LakeID)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
## -18434.2 -18381.1   9226.1 -18452.2     2693 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.6378 -0.2586 -0.0892  0.1320 21.0867 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 1.785e-02 0.133621
##  LakeID             (Intercept) 7.516e-03 0.086695
##  Sample_Month       (Intercept) 7.055e-04 0.026562
##  Sample_Year        (Intercept) 5.504e-04 0.023460
##  Residual                       6.858e-05 0.008282
## Number of obs: 2702, groups:  
## Sample_Year:LakeID, 349; LakeID, 69; Sample_Month, 6; Sample_Year, 6
## 
## Fixed effects:
##                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          -3.878588   0.120002 -32.321  < 2e-16 ***
## DreissenidsUninvaded -0.071907   0.080806  -0.890  0.37353    
## log(CA.SA)            0.018227   0.027491   0.663  0.50732    
## Mean_Depth_m         -0.019802   0.007044  -2.811  0.00493 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.640              
## log(CA.SA)  -0.464  0.072       
## Mean_Dpth_m -0.260 -0.009 -0.233
## convergence code: 0
## Model failed to converge with max|grad| = 0.433841 (tol = 0.001, component 1)
plot(fitted(Gfit), residuals(Gfit), xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2)
lines(smooth.spline(fitted(Gfit), residuals(Gfit)))

Model 2 - Gamma with Inverse Link

#model transformed
Gfit2<-glmer(TP_mg.L+.01 ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1 | Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df, family=Gamma(link="inverse"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.126138 (tol = 0.001, component 1)
summary(Gfit2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: TP_mg.L + 0.01 ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1 |  
##     Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 | Sample_Year:LakeID)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
## -21857.9 -21804.8  10938.0 -21875.9     2693 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7424 -0.4167 -0.1366  0.2090 14.5261 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  Sample_Year:LakeID (Intercept) 25.02369 5.0024  
##  LakeID             (Intercept) 15.42516 3.9275  
##  Sample_Month       (Intercept)  0.29972 0.5475  
##  Sample_Year        (Intercept)  0.69974 0.8365  
##  Residual                        0.05531 0.2352  
## Number of obs: 2702, groups:  
## Sample_Year:LakeID, 349; LakeID, 69; Sample_Month, 6; Sample_Year, 6
## 
## Fixed effects:
##                      Estimate Std. Error t value Pr(>|z|)    
## (Intercept)           44.9153     4.2107  10.667  < 2e-16 ***
## DreissenidsUninvaded   3.1243     3.0935   1.010 0.312522    
## log(CA.SA)            -0.7754     1.1621  -0.667 0.504594    
## Mean_Depth_m           0.9865     0.2959   3.334 0.000855 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.651              
## log(CA.SA)  -0.514  0.032       
## Mean_Dpth_m -0.258 -0.043 -0.253
## convergence code: 0
## Model failed to converge with max|grad| = 0.126138 (tol = 0.001, component 1)
plot(fitted(Gfit2), residuals(Gfit2), xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2)
lines(smooth.spline(fitted(Gfit2), residuals(Gfit2)))

Model 2b - Gamma with Inverse Link and lakes and year 2015 removed

#model transformed
Gfit2b<-glmer(TP_mg.L+.01 ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1|Sample_Year) + (1 | Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df4, family=Gamma(link="inverse"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0553478 (tol = 0.001, component 1)
summary(Gfit2b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: TP_mg.L + 0.01 ~ Dreissenids + log(CA.SA) + Mean_Depth_m + (1 |  
##     Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 | Sample_Year:LakeID)
##    Data: df4
## 
##      AIC      BIC   logLik deviance df.resid 
## -17857.3 -17806.3   8937.6 -17875.3     2105 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9642 -0.4499 -0.1464  0.2399 11.1878 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  Sample_Year:LakeID (Intercept) 14.76129 3.8420  
##  LakeID             (Intercept) 16.42993 4.0534  
##  Sample_Month       (Intercept)  0.17428 0.4175  
##  Sample_Year        (Intercept)  0.82420 0.9079  
##  Residual                        0.03906 0.1976  
## Number of obs: 2114, groups:  
## Sample_Year:LakeID, 272; LakeID, 65; Sample_Month, 6; Sample_Year, 5
## 
## Fixed effects:
##                      Estimate Std. Error t value Pr(>|z|)    
## (Intercept)           45.3284     4.7589   9.525  < 2e-16 ***
## DreissenidsUninvaded   2.5946     3.2172   0.806  0.41997    
## log(CA.SA)            -0.6284     1.2950  -0.485  0.62751    
## Mean_Depth_m           0.9417     0.3203   2.940  0.00328 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.634              
## log(CA.SA)  -0.544  0.058       
## Mean_Dpth_m -0.279 -0.025 -0.226
## convergence code: 0
## Model failed to converge with max|grad| = 0.0553478 (tol = 0.001, component 1)
plot(fitted(Gfit2b), residuals(Gfit2b), xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2)
lines(smooth.spline(fitted(Gfit2b), residuals(Gfit2b)))

Model 3 - Gamma with Log Link

-Errors: 1) Model is nearly unidentifiable 2) Model failed to converge -Tried re-scaling Mean Depth using log transformation

#model transformed
Gfit3<-glmer(TP_mg.L+.01 ~ Dreissenids + log(CA.SA) + log(Mean_Depth_m) + (1|Sample_Year) + (1 | Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df, family=Gamma(link="log"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.12433 (tol = 0.001, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
summary(Gfit3)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: TP_mg.L + 0.01 ~ Dreissenids + log(CA.SA) + log(Mean_Depth_m) +  
##     (1 | Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 |  
##     Sample_Year:LakeID)
##    Data: df
## 
##      AIC      BIC   logLik deviance df.resid 
## -21746.2 -21693.1  10882.1 -21764.2     2693 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.4888 -0.4107 -0.1374  0.1896 15.5550 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 0.0187517 0.13694 
##  LakeID             (Intercept) 0.0076974 0.08774 
##  Sample_Month       (Intercept) 0.0002307 0.01519 
##  Sample_Year        (Intercept) 0.0006536 0.02557 
##  Residual                       0.0565140 0.23773 
## Number of obs: 2702, groups:  
## Sample_Year:LakeID, 349; LakeID, 69; Sample_Month, 6; Sample_Year, 6
## 
## Fixed effects:
##                        Estimate Std. Error  t value Pr(>|z|)    
## (Intercept)          -3.5733901  0.0005725 -6241.66   <2e-16 ***
## DreissenidsUninvaded -0.0763196  0.0005726  -133.30   <2e-16 ***
## log(CA.SA)            0.0073759  0.0005727    12.88   <2e-16 ***
## log(Mean_Depth_m)    -0.1297226  0.0005727  -226.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv 0.001               
## log(CA.SA)  0.001  0.001        
## lg(Mn_Dpt_) 0.001  0.001  0.001 
## convergence code: 0
## Model failed to converge with max|grad| = 0.12433 (tol = 0.001, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
plot(fitted(Gfit3), residuals(Gfit3), xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2)
lines(smooth.spline(fitted(Gfit3), residuals(Gfit3)))

simGfit3  <-  simulateResiduals(Gfit3)
## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.
plotSimulatedResiduals(simGfit3) 
## plotSimulatedResiduals is deprecated, switch your code to using the plot function

Model 3b - Gamma with Log Link and lakes and Year 2015 removed

-Error: 1) Model failed to converge

#model transformed
Gfit3b<-glmer(TP_mg.L+.01 ~ Dreissenids + log(CA.SA) + log(Mean_Depth_m) + (1|Sample_Year) + (1 | Sample_Month) + (1|LakeID) + (1|Sample_Year:LakeID), data=df4, family=Gamma(link="log"))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00144156 (tol = 0.001, component 1)
summary(Gfit3b)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( log )
## Formula: TP_mg.L + 0.01 ~ Dreissenids + log(CA.SA) + log(Mean_Depth_m) +  
##     (1 | Sample_Year) + (1 | Sample_Month) + (1 | LakeID) + (1 |  
##     Sample_Year:LakeID)
##    Data: df4
## 
##      AIC      BIC   logLik deviance df.resid 
## -17841.4 -17790.5   8929.7 -17859.4     2105 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9626 -0.4515 -0.1518  0.2276 11.4646 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Sample_Year:LakeID (Intercept) 6.687e-03 0.081775
##  LakeID             (Intercept) 7.642e-03 0.087418
##  Sample_Month       (Intercept) 7.298e-05 0.008543
##  Sample_Year        (Intercept) 3.896e-04 0.019738
##  Residual                       3.907e-02 0.197663
## Number of obs: 2114, groups:  
## Sample_Year:LakeID, 272; LakeID, 65; Sample_Month, 6; Sample_Year, 5
## 
## Fixed effects:
##                       Estimate Std. Error t value Pr(>|z|)    
## (Intercept)          -3.659078   0.122282 -29.923   <2e-16 ***
## DreissenidsUninvaded -0.065173   0.070074  -0.930   0.3523    
## log(CA.SA)            0.007509   0.028380   0.265   0.7913    
## log(Mean_Depth_m)    -0.106227   0.041623  -2.552   0.0107 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) DrssnU l(CA.S
## DrssndsUnnv -0.583              
## log(CA.SA)  -0.536  0.052       
## lg(Mn_Dpt_) -0.568  0.076 -0.023
## convergence code: 0
## Model failed to converge with max|grad| = 0.00144156 (tol = 0.001, component 1)
plot(fitted(Gfit3b), residuals(Gfit3b), xlab = "Fitted Values", ylab = "Residuals")
abline(h = 0, lty = 2)
lines(smooth.spline(fitted(Gfit3b), residuals(Gfit3b)))

simGfit3b  <-  simulateResiduals(Gfit3b)
## Model family was recognized or set as continuous, but duplicate values were detected in the response. Consider if you are fitting an appropriate model.
plotSimulatedResiduals(simGfit3b) 
## plotSimulatedResiduals is deprecated, switch your code to using the plot function